#include "Gas.h"

int indx(int i, int j){
  return (i*DIM+j);
};

Gas::Gas(int const N){
	this->N= N;
  this->theta = new double[N];
  this->omega = new double[N];
  
  this->vec_m = new double[N*DIM];
  this->vec_M = new double[DIM];
  this->force = new double[N];
  
  this->Ek =0.0;
  this->Ep =0.0;
}

Gas::~Gas(){
  delete[] theta;
  delete[] omega;
  delete[] vec_M;
  delete[] vec_m;
  delete[] force;
};

void Gas::initCoord(){
  double pi=3.14159;
  for (int i=0; i<this->N ; i++){
    this->theta[i]=(2.0*(i+0.5)/(this->N) -1 )*pi ;
  }
};

void Gas::initVeloc(){
  for (int i=0; i<this->N ; i++){
    this->omega[i]=0.0;
  }
};


// Calculate functions
void Gas::calculateVec_m(){
  for (int i=0; i< this->N; i++){
    this->vec_m[indx(i,0)] = cos(this->theta[i]);
    this->vec_m[indx(i,1)] = sin(this->theta[i]);
  }
}

void Gas::calculateVec_M(){
  this->vec_M[0]=0;
  this->vec_M[1]=0;
  for (int i=0; i<this->N; i++){
    this->vec_M[0] += this->vec_m[indx(i,0)];
    this->vec_M[1] += this->vec_m[indx(i,1)];
  }

  this->vec_M[0] = this->vec_M[0]/this->N;
  this->vec_M[1] = this->vec_M[1]/this->N;
//  cout<<this->vec_M[0]<<" "<<this->vec_M[1]<<endl;
}

void Gas::calculateForce(double eps){
  for (int i=0; i<this->N; i++){
    this->force[i] = -eps*(this->vec_M[0]*sin(this->theta[i]) - this->vec_M[1]*cos(this->theta[i]));
  }
}

// Energy
void Gas::calculateEkin(){
	this->Ek=0.0;
	for (int i=0; i<this->N; i++){
		this->Ek += this->omega[i]*this->omega[i];
	}
	this->Ek = 0.5*this->Ek;
}

void Gas::calculateEpot(double eps){
	this->Ep=0.0;
	for (int i=0; i<this->N; i++){
		for (int j=0; j<this->N; j++){
			this->Ep += 1
				-(
					 this->vec_m[indx(i,0)]*this->vec_m[indx(j,0)]
					+this->vec_m[indx(i,1)]*this->vec_m[indx(j,1)]
				 );
		}
	}
	this->Ep = (eps/(2.0*N))*this->Ep;
}

// Update functions
void Gas::updateCoord(Gas* oldgas, double dt, double dt2){
	for (int i=0; i<this->N; i++){
    this->theta[i] = oldgas->theta[i]+oldgas->omega[i]*dt+0.5*dt2*oldgas->force[i];
  }
}

void Gas::updateVeloc(Gas* oldgas, double dt){
  for (int i=0; i<this->N ; i++){
    this->omega[i]=oldgas->omega[i]+ 0.5*dt*(this->force[i]+oldgas->force[i]);
  }
}

// Save data to file
void Gas::writeMacroState(ofstream *ptrfile, int ttime, double dt){
	(*ptrfile) << dt*ttime<<" ";
	(*ptrfile) << this->Ek << " ";
	(*ptrfile) << this->Ep << " ";
	(*ptrfile) << this->vec_M[0] << " ";
	(*ptrfile) << this->vec_M[1] << " ";
	(*ptrfile) << endl;
}

void Gas::writeMicroState(ofstream *ptrfile, int ttime, double dt){
	(*ptrfile) << dt*ttime<<" ";
	for (int i =0; i<this->N; i++){
		(*ptrfile) << this->omega[i] << " ";
		(*ptrfile) << this->theta[i] << " ";	
	}
	(*ptrfile) << endl;
}

// print functions
void Gas::printCoord(){
	for (int i=0; i<this->N ; i++){
    cout<<this->theta[i]<<" ";
  }
  cout<<endl;
}

void Gas::printVeloc(){
  for (int i=0; i<this->N ; i++){
    cout << this->omega[i]<<" ";
  }
  cout<<endl;
}

void Gas::printVec_m(){
  cout<<"_____ Begin printVec_m() "<<endl;
  for (int i=0; i< this->N; i++){
  	cout<< this->vec_m[indx(i,0)]<< ", " << this->vec_m[indx(i,1)] <<endl;
  }
  cout<<"_____ End printVec_m() "<<endl;
}

void Gas::printVec_M(){
  cout<<"_____ Begin printVec_M() "<<endl;
	cout<< this->vec_M[0]<< ", " << this->vec_M[1] <<endl;
  cout<<"_____ End printVec_M() "<<endl;
}

void Gas::printForce(){
	for (int i=0; i<this->N; i++){
		cout<< this->force[i] <<endl;
	}
}

